//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
/*************************************************************************
 * template<class Mat_t>
 * class  LapackSymAllSolver: public EigenSolverBase<Mat_t>
 * using a simple format for the storage
 *
 * solving ALL eigen value problems using LAPACK routine
 *************************************************************************/
#ifndef OHMMS_LAPACKSYMALLSOLVERALL_H
#define OHMMS_LAPACKSYMALLSOLVERALL_H

#ifndef EIGENSPACEDIMENSION_H
#include "EigSolver/EigenSpaceDimension.h"
#endif
#include "OhmmsNumeric/OhmmsBlasf.h"

#ifdef ADD_
#define dsyev dsyev_
#define zheev zheev_
#define dsygv dsygv_
#define zhegv zhegv_
#endif

extern "C" {
  void dsyev(const char&, const char&, const int&,

             double*, const int&, double*, double*, int&, int&);

  void zheev(const char&, const char&, const int&,
             complex<double>* a, const int& lda, double* w, 
             complex<double>* work, int& lwork, double* rwork, int& info);

  void dsygv(const int&, const char&, const char&, const int& , 
             double* , const int& , double* , const int& , double* , 
             double* , int& , int&);
  void zhegv(const int&, const char&, const char&, 
             const int& , complex<double>* , const int& , complex<double>* , 
             const int& , double* , 
             complex<double>* , int& , double*, int&);

}
template<class Mat_t, class EigVec_t, class EigVal_t>
class LapackSymSolver {

  EigenSpaceDimension& States;       

public:

  //typedef typename Mat_t::Element_t     Element_t;
  typedef typename Mat_t::Type_t        Element_t;
  typedef typename EigVal_t::value_type real_t;
  typedef vector<Element_t>             WorkArray_t;
  typedef vector<real_t>                RealArray_t;

public:

  LapackSymSolver(EigenSpaceDimension& dim): 
    States(dim), N(0), Jobz('V'), Uplo('L'), Itype(1) { }

  ~LapackSymSolver() { }

  bool getEigvectors(Mat_t& H, EigVec_t&, EigVal_t& );
  bool getEigvalues(Mat_t& H, EigVec_t&, EigVal_t& );

  bool getEigvectors(Mat_t& H, Mat_t& S, EigVec_t&, EigVal_t& );
  bool getEigvalues(Mat_t& H, Mat_t& S,  EigVec_t&, EigVal_t& );

  int create(Mat_t& H, Mat_t& S, Mat_t& Wfs) {
    H.resize(States.N, States.N);
    S.resize(States.N, States.N);
    return 0;
  }

  int create(Mat_t& H, Mat_t& Wfs) {
    H.resize(States.N, States.N);
    return 0;
  }

private:
  int N;
  int lwork;
  int info;
  int Itype;
  char Jobz, Uplo;
  WorkArray_t workd; 
  RealArray_t rworkd;
  RealArray_t eigVal;

  void initializeWorkSpace();

  void solve(double* a) {
    //dsyev(Jobz, Uplo, States.N, a, States.N, eigVal.begin(), 
    //      workd.begin(), lwork, info);
    dsyev(Jobz, Uplo, States.N, a, States.N, &(eigVal[0]), 
          &(workd[0]), lwork, info);
  }

  void solve(double* a, double *b) {
    dsygv(Itype, Jobz, Uplo, States.N, a, States.N, b, 
          States.N, &(eigVal[0]), &(workd[0]), lwork, info);
  }

  void solve(complex<double>* a) {
    if(rworkd.empty()) { rworkd = RealArray_t(3*N); }
    zheev(Jobz, Uplo, States.N, a, States.N, &(eigVal[0]), 
          &(workd[0]), lwork, &(rworkd[0]), info);
  }


  void solve(complex<double>* a, complex<double>* b) {
    if(rworkd.empty()) { rworkd = RealArray_t(3*N); }
    zhegv(Itype, Jobz, Uplo, States.N, a, States.N, b, States.N, 
          &eigVal[0], &workd[0], lwork, &rworkd[0],info);

  }

};

template<class Mat_t, class EigVec_t, class EigVal_t>
bool 
LapackSymSolver<Mat_t,EigVec_t,EigVal_t>::getEigvectors(Mat_t& H, 
                                         EigVec_t& eigvec, 
                                         EigVal_t& eigval) {
  initializeWorkSpace();

  Jobz = 'V';
  solve(H.begin()); 
  if(info < 0) return false;
  else for(int i=0; i<States.M; i++) eigval[i] = eigVal[i];
  return true;
}

template<class Mat_t, class EigVec_t, class EigVal_t>
bool 
LapackSymSolver<Mat_t,EigVec_t,EigVal_t>::getEigvalues(Mat_t& H,
                                                       EigVec_t& eigvec, 
                                                       EigVal_t& eigval) {
  initializeWorkSpace();
  Jobz = 'N';
  solve(H.begin());    

  if(info < 0) return false;
  else for(int i=0; i<States.M; i++) eigval[i] = eigVal[i];
  return true;
}

template<class Mat_t, class EigVec_t, class EigVal_t>
bool 
LapackSymSolver<Mat_t,EigVec_t,EigVal_t>::getEigvectors(Mat_t& H, Mat_t& S, 
                                                        EigVec_t& eigvec, 
                                                        EigVal_t& eigval) {
  initializeWorkSpace();
  Jobz = 'V';
  solve(H.begin(),S.begin());            // solve eigen vector problems
  if(info < 0) {
    cout << "Error has occurred. info = " << info << endl; 
    return false;
  }
  else {
    for(int i=0; i<States.M; i++) eigval[i] = eigVal[i];
  }
  return true;
}

template<class Mat_t, class EigVec_t, class EigVal_t>
bool 
LapackSymSolver<Mat_t,EigVec_t,EigVal_t>::getEigvalues(Mat_t& H, Mat_t& S, 
                                                       EigVec_t& eigvec, 
                                                       EigVal_t& eigval) {
  initializeWorkSpace();
  Jobz = 'N';
  solve(H.begin(), S.begin());          // solve eigen value problems

  if(info < 0) return false;
  else for(int i=0; i<States.M; i++) eigval[i] = eigVal[i];
  return true;
}

template<class Mat_t, class EigVec_t, class EigVal_t>
void
LapackSymSolver<Mat_t,EigVec_t,EigVal_t>::initializeWorkSpace() {

  if(N!= States.N) {
    N = States.N;
    lwork = 3*N;
    if(workd.size() < lwork) {
      workd.erase(workd.begin(), workd.end());
      workd = WorkArray_t(lwork);
    }
    if(eigVal.size() < N) {
      eigVal.erase(eigVal.begin(), eigVal.end());
      eigVal = RealArray_t(N);
    }
  }
}
#endif

/***************************************************************************
 * $RCSfile: LapackSymSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: LapackSymSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

